home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-ppc folder / startup / fft.i < prev    next >
Text File  |  1995-07-26  |  11KB  |  320 lines
  1. /*
  2.    FFT.I
  3.    Yorick interface to Swarztrauber FFT routines, modified for strides.
  4.  
  5.    $Id: fft.i,v 1.1 1993/08/06 23:13:12 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* ------------------------------------------------------------------------ */
  11.  
  12. func fft(x, ljdir, rjdir, setup=)
  13. /* DOCUMENT fft(x, direction)
  14.             fft(x, ljdir, rjdir)
  15.          or fft(x, ljdir, rjdir, setup=workspace)
  16.      returns the complex Fast Fourier Transform of array X.
  17.      The DIRECTION determines which direction the transform is in --
  18.      e.g.- from time to frequency or vice-versa -- as follows:
  19.  
  20.      DIRECTION    meaning
  21.      ---------    -------
  22.          1        "forward" transform (coefficients of exp(+i * 2*pi*kl/N))
  23.               on every dimension of X
  24.         -1        "backward" transform (coefficients of exp(-i * 2*pi*kl/N))
  25.               on every dimension of X
  26.      [1,-1,1]     forward transform on first and third dimensions of X,
  27.                   backward transform on second dimension of X (any other
  28.           dimensions remain untransformed)
  29.      [-1,0,0,1]   backward transform on first dimension of X, forward
  30.                   transform on fourth dimension of X
  31.     etc.
  32.  
  33.      The third positional argument, if present, allows the direction
  34.      of dimensions of X to be specified relative to the final dimension
  35.      of X, instead of relative to the first dimension of X.  In this
  36.      case, both LJDIR and RJDIR must be vectors of integers -- the
  37.      scalar form is illegal:
  38.  
  39.         LJDIR    RJDIR      meaning
  40.         -----    -----      -------
  41.     []        [1]       forward transform last dimension of X
  42.     [1]        []       forward transform first dimension of X
  43.     []        [-1,-1]   backward transform last two dimensions of X,
  44.                         leaving any other dimensions untransformed
  45.      [-1,0,0,1]    []       backward transform on first dimension of X,
  46.                             forward transform on fourth dimension of X
  47.         []      [-1,0,0,1]  backward transform on 4th to last dimension of X,
  48.                             forward transform on last dimension of X
  49.     etc.
  50.  
  51.      Note that the final element of RJDIR corresponds to the last dimension
  52.      of X, while the initial element of LJDIR corresponds to the first
  53.      dimension of X.
  54.  
  55.      The explicit meaning of "forward" transform -- the coefficients of
  56.      exp(+i * 2*pi*kl/N) -- is:
  57.  
  58.      result for j=1,...,n
  59.  
  60.                 result(j)=the sum from k=1,...,n of
  61.  
  62.                       x(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
  63.  
  64.                             where i=sqrt(-1)
  65.  
  66.      Note that the result is unnormalized.  Applying the "backward"
  67.      transform to the result of a "forward" transform returns N times
  68.      the original vector of length N.  Equivalently, applying either
  69.      the "forward" or "backward" transform four times in succession
  70.      yields N^2 times the original vector of length N.
  71.  
  72.      Performing the transform requires some WORKSPACE, which can be
  73.      set up beforehand by calling fft_setup, if fft is to be called
  74.      more than once with arrays X of the same shape.  If no setup
  75.      keyword argument is supplied, the workspace allocation and setup
  76.      must be repeated for each call.
  77.  
  78.    SEE ALSO: roll, fft_setup, fft_inplace
  79.  */
  80. {
  81.   result= complex(x);
  82.   fft_inplace, result, ljdir, rjdir, setup=setup;
  83.   return result;
  84. }
  85.  
  86. func fft_inplace(x, ljdir, rjdir, setup=)
  87. /* DOCUMENT fft_inplace, x, direction
  88.          or fft_inplace, x, ljdir, rjdir
  89.      or fft_inplace, x, ljdir, rjdir, setup=workspace
  90.      is the same as the fft function, except that the transform is
  91.      performed "in_place" on the array X, which must be of type complex.
  92.    SEE ALSO: fft, fft_setup
  93.  */
  94. {
  95.   if (structof(x)!=complex) error, "expecting complex argument";
  96.  
  97.   /* form list of dimension lengths (dims) and corresponding list of
  98.      transform directions (dirs) */
  99.   dims= dimsof(x);
  100.   ndims= dims(1);
  101.   if (ndims<1) return;
  102.   dirs= fft_dirs(ndims, ljdir, rjdir);
  103.   list= where(dirs);  /* list of dimensions to be transformed */
  104.   n= numberof(list);
  105.   if (!n) return;
  106.  
  107.   /* setup the workspace now if it hasn't already been done */
  108.   if (is_void(setup)) setup= fft_setup(dims, dirs);
  109.  
  110.   /* get the strides associated with each dimension */
  111.   dims= dims(2:0);
  112.   stds= array(1, ndims);
  113.   for (i=1 ; i<ndims ; i++) stds(i+1)= stds(i)*dims(i);
  114.   tops= numberof(x)/(stds*dims);
  115.  
  116.   /* do the requested transforms */
  117.   nlist= *setup(1);
  118.   setup= *setup(2);
  119.   if (is_void(setup)) error, "no workspace for any dimension length";
  120.   for (i=1 ; i<=n ; i++) {
  121.     j= list(i);
  122.     ndims= dims(j);  /* clobbers original ndims */
  123.     ws= setup(where(nlist==ndims)(1));
  124.     if (is_void(ws)) error, "no workspace for dimension length "+pr1(ndims);
  125.     fft_raw, dirs(j), x, stds(j), ndims, tops(j), ws;
  126.   }
  127. }
  128.  
  129. func fft_setup(dims, ljdir, rjdir)
  130. /* DOCUMENT workspace= fft_setup(dimsof(x))
  131.          or workspace= fft_setup(dimsof(x), direction)
  132.          or workspace= fft_setup(dimsof(x), ljdir, rjdir)
  133.      allocates and sets up the workspace for a subsequent call to
  134.             fft(X, DIRECTION, setup=WORKSPACE)
  135.      or
  136.             fft(X, LJDIR, RJDIR, setup=WORKSPACE)
  137.      The DIRECTION or LJDIR, RJDIR arguments compute WORKSPACE only for
  138.      the dimensions which will actually be transformed.  If only the
  139.      dimsof(x) argument is supplied, then WORKSPACE will be enough to
  140.      transform any or all dimensions of X.  With DIRECTION or LJDIR, RJDIR
  141.      supplied, WORKSPACE will only be enough to compute the dimensions
  142.      which are actually to be transformed.  The WORKSPACE does not
  143.      depend on the sign of any element in the DIRECTION (or LJDIR, RJDIR),
  144.      so you can use the same WORKSPACE for both "forward" and "backward"
  145.      transforms.
  146.  
  147.      Furthermore, as long as the length of any dimensions of the array
  148.      X to be transformed are present in WORKSPACE, it may be used in
  149.      a call to fft with the array.  Thus, if X were a 25-by-64 array,
  150.      and Y were a 64-vector, the following sequence is legal:
  151.           ws= fft_setup(dimsof(x));
  152.       xf= fft(x, 1, setup=ws);
  153.       yf= fft(y, -1, setup=ws);
  154.  
  155.      The WORKSPACE required for a dimension of length N is 6*N+15 doubles.
  156.  
  157.    SEE ALSO: fft, dimsof, fft_inplace
  158.  */
  159. {
  160.   ndims= dims(1);
  161.   if (ndims<1) return array(pointer, 2);
  162.   dirs= fft_dirs(ndims, ljdir, rjdir);
  163.   list= where(dirs);  /* list of dimensions to be transformed */
  164.   n= numberof(list);
  165.   if (!n) return array(pointer, 2);
  166.   dims= dims(1+list);
  167.   if (n>1) {
  168.     dims= dims(sort(dims));
  169.     list= grow([0],dims)(dif);
  170.     /* eliminate duplicate dimension lengths */
  171.     dims= dims(where(list));
  172.     n= numberof(dims);
  173.   }
  174.   setup= array(pointer, n);
  175.   for (i=1 ; i<=n ; i++) {
  176.     ndims= dims(i);
  177.     ws= array(0.0, 6*ndims+15);
  178.     fft_init, ndims, ws;
  179.     setup(i)= &ws;
  180.   }
  181.   return [&dims, &setup];
  182. }
  183.  
  184. func fft_dirs(ndims, ljdir, rjdir)
  185. {
  186.   nolj= is_void(ljdir);
  187.   norj= is_void(rjdir);
  188.   if (nolj && norj) return array(1, ndims);  /* mostly for fft_setup --
  189.                             also makes fft(x) work */
  190.   if (nolj || dimsof(ljdir)(1)) {
  191.     /* detailed direction list present */
  192.     dirs= array(0, ndims);          /* all 0 initially */
  193.     if (!nolj) {                    /* fill in left-justified directions */
  194.       list= where(ljdir);
  195.       dirs(list)= ljdir(list);
  196.     }
  197.     if (!norj) {                    /* fill in right-justified directions */
  198.       list= where(rjdir);
  199.       dirs(list-numberof(rjdir) + ndims)= rjdir(list);
  200.     }
  201.   } else {
  202.     /* all dimensions get the same direction */
  203.     dirs= array(ljdir, ndims);
  204.   }
  205.   return dirs;
  206. }
  207.  
  208. extern fft_init;
  209. /* PROTOTYPE
  210.    void cffti(long n, double array ws)
  211.  */
  212. /* DOCUMENT fft_init, n, wsave
  213.      Swarztrauber's cffti.  This actually requires wsave=array(0.0, 4*n+15),
  214.      instead of the 6*n+15 doubles of storage used by fft_raw to handle the
  215.      possibility of multidimensional arrays.  If the storage matters, you
  216.      can call cfftf and/or cfftb as the Yorick functions fft_fraw and/or
  217.      fft_braw.
  218.  */
  219.  
  220. extern fft_fraw;
  221. /* PROTOTYPE
  222.    void cfftf(long n, complex array c, double array ws)
  223.  */
  224. /* DOCUMENT fft_fraw, n, c, wsave
  225.      Swarztrauber's cfftf.  You can use this to avoid the additional
  226.      2*N storage incurred by fft_setup.
  227.  */
  228.  
  229. extern fft_braw;
  230. /* PROTOTYPE
  231.    void cfftb(long n, complex array c, double array ws)
  232.  */
  233. /* DOCUMENT fft_braw, n, c, wsave
  234.      Swarztrauber's cfftb.  You can use this to avoid the additional
  235.      2*N storage incurred by fft_setup.
  236.  */
  237.  
  238. extern fft_raw;
  239. /* PROTOTYPE
  240.    void cfft2(long idir, complex array c, long istd, long n, long n2,
  241.               pointer ws)
  242.  */
  243.  
  244. /* ------------------------------------------------------------------------ */
  245.  
  246. func roll(x, ljoff, rjoff)
  247. /* DOCUMENT roll(x, ljoff, rjoff)
  248.          or roll, x, ljoff, rjoff
  249.      or roll(x)
  250.      or roll, x
  251.  
  252.      "rolls" selected dimensions of the array X.  The roll offsets
  253.      LJOFF and RJOFF (both optional) work in the same fashion as the
  254.      LJDIR and RJDIR arguments to the fft function:
  255.  
  256.         A scalar LJDIR (and nil RJDIR) rolls all dimensions of X by
  257.         the specified offset.
  258.     Otherwise, the elements of the LJDIR vector [ljoff1, ljoff2, ...]
  259.     are used as the roll offsets for the first, second, etc.
  260.     dimensions of X.
  261.     Similarly, the elements of the RJDIR vector [..., rjoff1, rjoff0]
  262.     are matched to the final dimensions of X, so the next to last
  263.     dimension is rolled by rjoff1 and the last dimension by rjoff0.
  264.  
  265.     As a special case (mostly for use with the fft function), if
  266.     both LJDIR and RJDIR are nil, every dimension is rolled by
  267.     half of its length.  Thus,
  268.        roll(x)
  269.     it equivalent to
  270.        roll(x, dimsof(x)(2:0)/2)
  271.  
  272.      The result of the roll function is complex if X is complex, otherwise
  273.      double (i.e.- any other array type is promoted to type double).  If
  274.      roll is invoked as a subroutine, the operation is performed in place.
  275.  
  276.    SEE ALSO: fft
  277.  */
  278. {
  279.   /* get array to be transformed */
  280.   if (structof(x)==complex) {
  281.     stds1= 2;
  282.     if (!am_subroutine()) x= complex(x);  /* copy to avoid clobbering */
  283.   } else {
  284.     stds1= 1;
  285.     if (!am_subroutine()) x= double(x);  /* copy to avoid clobbering */
  286.     else if (structof(x)!=double)
  287.       error, "in-place roll needs double or complex argument";
  288.   }
  289.   dims= dimsof(x);
  290.   ndims= dims(1);
  291.   n= numberof(x);
  292.   if (ndims<1 || n<2) return x;
  293.   dims= dims(2:0);
  294.  
  295.   /* get offset arguments */
  296.   if (is_void(ljoff) && is_void(rjoff)) ljoff= dims/2;
  297.   else ljoff= fft_dirs(ndims, ljoff, rjoff);
  298.   ljoff%= dims;
  299.   if (noneof(ljoff)) return x;
  300.  
  301.   /* get strides */
  302.   stds= array(stds1, ndims);   /* note extra factor of 2 for complex */
  303.   for (i=1 ; i<ndims ; i++) stds(i+1)= stds(i)*dims(i);
  304.   tops= (stds1*n)/(stds*dims);
  305.  
  306.   /* perform roll operation */
  307.   px= &x;   /* avoids complex<-->double type conversion problems */
  308.   ws= array(0.0, max(dims(where(ljoff))));
  309.   for (i=1 ; i<=ndims ; i++)
  310.     _roll2, px, ljoff(i), stds(i), dims(i), tops(i), ws;
  311.  
  312.   return x;
  313. }
  314.  
  315. extern _roll2;
  316. /* PROTOTYPE
  317.    void roll2(pointer a, long ioff, long istd, long n, long n2,
  318.               double array ws)
  319.  */
  320.